#include "coccurance.h"

// Param list contains height, width, gray scale level and
// distance between two neighbors. Default value for distance is 1 which
// stands for direct neighbor of each pixel
Coccurance::Coccurance(int h, int w, int level) :
    height(h), width(w), nGray(level)
{
    delta1 = delta2 = delta3 = 0.0;

    init();
}

Coccurance::~Coccurance()
{
    for(int i = 0; i < 4; i++)
    {
        for(int p = 0; p < nGray; p++)
        {
            for(int q = 0; q <= p; q++)
            {
                delete glchs[i][p][q].listnode;
            }
            delete[] glchs[i][p];
        }
        delete[] glchs[i];
        delete head[i];
    }

    delete[] glchs;
    delete[] head;
    delete[] tail;
    delete stats;
}

/*****************************
    HashNode*** glchs;
    ListNode** head;
    ListNode** tail;
    Stat* stats;

    Since symmetrical GLCM is considered, only to initialize lower half
    of 'glchs' matrix. Also one 'glchs' for each direction. Dimension of
    each 'glchs' is based on number of gray level 'nGray'. Each cell of
    'glchs' stores a listnode that contains (i, j) as well as coccurrance
    prob of(i, j)
*******************************/
void Coccurance::init()
{
    glchs = new HashNode** [4];
    head = new ListNode* [4];
    tail = new ListNode* [4];
    stats = new Stat();

    for(int m = 0; m < 4; m++)
    {
        glchs[m] = new HashNode* [nGray];

        head[m] = new ListNode();
        tail[m] = head[m];

        // Initialize 2-D hash table
        for(int i = 0; i < nGray; i++)
        {
            // Each row has increasing capacities to store coccurance
            glchs[m][i] = new HashNode[i+1];

            // Symetrical matrix is used here
            // Only lower half of matrix is used to store nodes'
            // probabilities
            for(int j = 0; j <= i; j++)
            {
                ListNode* node = new ListNode();
                node->x1 = i;
                node->x2 = j;
                node->next = 0;
                glchs[m][i][j].prob = 0.0;
                glchs[m][i][j].listnode = node;
            }
        }

        head[m]->next = glchs[m][0][0].listnode;
    }   
}

//  Reset all variables, so each cell within 'glchs' will not
//  point to any other nodes. Then reset head and tail of all
// four directions.
void Coccurance::reset()
{    
    for(int m = 0; m < 4; m++)
    {
        for(int i = 0; i < nGray; i++)
        {
            // Symmetrical matrix is used here
            for(int j = 0; j <= i; j++)
            {
                glchs[m][i][j].listnode->next = 0;
                glchs[m][i][j].prob = 0.0;
            }
        }

        head[m]->next = 0;
        tail[m] = head[m];
    }
    stats->con = 0.0;
    stats->dis = 0.0;
    stats->uni = 0.0;
    stats->idm = 0.0;
    stats->ent = 0.0;
    stats->cor = 0.0;
    stats->inv = 0.0;
    stats->var = 0.0;
}


// Passing address of first element pointed by 'data'
// as well as distance 'd' between 2 neighbor pixels
void Coccurance::setParam(int *data, int d)
{
    // Prerequisites for calculation of coccurrance
    if(data == 0 || d >= height || d >= width || d <= 0) return;

    // Total number of coccuring pairs
    // Symmetrical matrics is considered
    // delta1 is for 0
    // delta2 is for 45 and 135
    // delta3 is for 90
    delta1 = 1.0 / ((float) height*(width-d));
    delta2 = 1.0 / ((float) (width-d)*(height-d));
    delta3 = 1.0 / ((float) (height-d)*width);

    // build coccurance matrix for int *d
    build_list(data, d);
    // calculate texture features
    calc_stat();
}

// Build coccurrance for two neighbor pixels at distance 'd'
// There are four directons to be considered neighbor:
// 0, 45, 90 and 135 degree. 'data' is enclosed in region of
// size (height x width)
//
void Coccurance::build_list(int *data, int d)
{
    // i - height
    // j - width
    int i, j;

    for(i = d; i < height; i++)
    {
        for(j = d; j < width-d; j++)
        {
            find_insert(*(data+i*width+j), *(data+i*width+j+d), 0);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j+d), 45);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j), 90);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j-d), 135);
        }
    }

    for(j = 0; j < d; j++)
    {
        for(i = d; i < height; i++)
        {
            find_insert(*(data+i*width+j), *(data+i*width+j+d), 0);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j+d), 45);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j), 90);
        }
    }

    for(i = 0; i < d; i++)
    {
        for(j = 0; j < width-d; j++)
        {
            find_insert(*(data+i*width+j), *(data+i*width+j+d), 0);
        }
    }

    for(j = width-d; j < width; j++)
    {
        for(i = d; i < height; i++)
        {
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j), 90);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j-d), 135);
        }
    }
}

/*
 *  Given a pair (x1, x2) at 'angle', find out the probabilities
 *  between them or cooocurance. If it is 0, then create a new
 *  node and append to the end of linked list. If it is not 0,
 *  simply update the probabilities of existed node in glchs matrix.
 *
 */
void Coccurance::find_insert(int x1, int x2, int angle)
{
    int index = angle / 45;

    // Quantize value down to 'ng' level
    x1 = (int) x1*nGray/256;
    x2 = (int) x2*nGray/256;

    // Force x1 >= x2 for lower half matrix
    if(x1 < x2)
    {
        int temp = x1;
        x1 = x2;
        x2 = temp;
    }

    float prob = glchs[index][x1][x2].prob;
    ListNode* n = glchs[index][x1][x2].listnode;

    // If prob==0.0, then append it to the end of sequence
    // pointed by tail pointer. This is because that this
    // pair of pixels/elements has not been encountered
    if(prob == 0.0)
    {
        tail[index]->next = n;
        tail[index] = n;
    }

    if(angle == 0)
    {
        glchs[index][x1][x2].prob = prob + delta1;
    }
    else if(angle == 45)
    {
        glchs[index][x1][x2].prob = prob + delta2;
    }
    else if(angle == 90)
    {
        glchs[index][x1][x2].prob = prob + delta3;
    }
    else if(angle == 135)
    {
        glchs[index][x1][x2].prob = prob + delta2;
    }
}


// Calculate statistics for current WSIZE x WSIZE
// block. Starting from the head node. Final value
// will be average of all four directions.
void Coccurance::calc_stat()
{
    ListNode* t;

    // Contrast, Uniformity, Dissimilarity, Entropy
    // Inverse Difference Moment, Correlation,
    // Inverse Difference and covariance
    float con=0.0, uni=0.0, dis=0.0, idm=0.0;
    float ent=0.0, cor=0.0, inv=0.0, var=0.0;
    float avg=0.0, variance=0.0, prob=0.0;

    int x1 = 0, x2 = 0;

    // Go through each angle, 0, 45, 90 and 135
    for(int i = 0; i < 4; i++)
    {
        // Obtain first listnode of each direction
        t = head[i]->next;

        // Reset average and covariance values
        avg = variance = 0.0;

        // Get average and covariance for current linked list
        if(t != 0)
        {
            getAvgVar(i, t, avg, variance);
            var += variance;
        }

        // Iterate through the linked list
        while(t != 0)
        {
            x1 = t->x1;
            x2 = t->x2;
            prob = glchs[i][x1][x2].prob;

            con += prob * (x1-x2)*(x1-x2);
            dis += prob * abs(float (x1-x2));
            uni += prob * prob;
            idm += prob / (1+(x1-x2)*(x1-x2));
            ent += prob * (log(prob));
            cor += (x1-avg)*(x2-avg) * prob/variance;
            inv += prob / (1+x1-x2);

            t = t->next;
        }
    }

    // Obtain average across all four directions
    stats->con = con / 4.0;
    stats->uni = uni / 4.0;
    stats->dis = dis / 4.0;
    stats->idm = idm / 4.0;
    stats->ent = -1*ent / 4.0;
    stats->cor = cor / 4.0;
    stats->inv = inv / 4.0;
    stats->var = var / 4.0;

 }


// Return calculated feature value based
// on 'func' which is 0-based.
float Coccurance::avg_stat(int func)
{
    float value = 0.0;

    switch (func)
    {
    case 0 :    // Contrast
        value = stats->con;
        break;
    case 1 :    // Dissimilarity
        value = stats->dis;
        break;
    case 2 :    // Uniformity
        value = stats->uni;
        break;
    case 3 :    // Inverse difference moment
        value = stats->idm;
        break;
    case 4 :    // Entropy
        value = stats->ent;
        break;
    case 5 :    // Correlation
        value = stats->cor;
        break;
    case 6 :    // Inverse difference
        value = stats->inv;
        break;
    case 7 :    // Covariance
        value = stats->var;
        break;
    default:
        break;

    }

    return value;
}

/*
*   Get average and covariance for linked list. The node
*   is the head node of each direction indicated by index
*
*/
void Coccurance::getAvgVar(int index, ListNode *node, float &avg, float &var)
{
    int x1 = 0,  x2 = 0;
    float prob = 0.0;
    ListNode *t = node;

    // Calculate average of the linked list
    while(t != 0)
    {
        x1 = t->x1;
        x2 = t->x2;
        prob = glchs[index][x1][x2].prob;

        avg += x1 * prob;

        t = t->next;
    }

    // Reset t to  point to start of linked list
    t = node;

    // Calculate covariance of the linked list
    while(t != 0)
    {
        x1 = t->x1;
        x2 = t->x2;
        prob = glchs[index][x1][x2].prob;

        var += (x1 - avg) * (x1 - avg) * prob;

        t = t->next;
    }

    t = 0;
}
